function [x,z,return_data,transect_dist] = rslice_xslice (vslice_x_pt, vslice_y_pt, grid, gdata )
% RSLICE_XSLICE:  extract the horizontal and vertical components of the vertical slice
%
% PARAMETERS:
% Input:
%    vslice_x_pt, vslice_y_pt:
%        Arrays of x and y points within the grid.
%        system.
%     grid:
%        Structure with four fields, x, y, z, and mask.
%     gdata:
%        
% Output:
%     x:
%        2D horizontal dimension component of the slice.  
%     z:
%        2D vertical dimension component of the slice (depth)
%     return_data:
%        2D data component of the slice.  
%     transect_dist:
%        1D horizontal distance along the slice.
%
% Calls:
%     rslice_retrieve_data
% Called By:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% $Name: rslice-0_9_2 $
% $Id: rslice_xslice.m,v 1.8 2006/05/03 21:27:34 jevans Exp $
% AUTHOR:  johnevans@acm.org
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


x = [];
y = [];
u = [];
transect_dist = [];

if isempty ( gdata.z_type )
	z_type = '2D';
else
	z_type = '3D';
end


%
% We need to find out if any of the slice points are outside of
% the grid.  Data for these points shall be NaN.
% The interpolation scheme may think otherwise, but we know better.
[r,c] = size(grid.x);
boundary_x = [grid.x(1:r-1,1); grid.x(r,1:c-1)'; grid.x(r:-1:2,c); grid.x(1,c:-1:2)'];
boundary_y = [grid.y(1:r-1,1); grid.y(r,1:c-1)'; grid.y(r:-1:2,c); grid.y(1,c:-1:2)'];
vslice_pt_inside_grid = inpolygon ( vslice_x_pt, vslice_y_pt, boundary_x, boundary_y );


%
% Where is the mask equal to zero?  Use this information to NaN out
% land spots down the water column.
land_inds = find(grid.mask == 0);
grid.mask(land_inds) = NaN;


%
% retrieve the data at the current timestep
data = rslice_retrieve_data ( gdata );

%
% If there are more than two dimension, then z should be the leading dimension.
if ( length(size(data)) > 2 )
	nz = size(data,1);
else
	nz = 1;
end


if ~isempty ( grid.mask )
	%
	% Nan out land indices
	if isempty ( gdata.z_type )
		mask = grid.mask;
	else
		%
		% ZYX data has to be masked at all z levels
		mask = repmat ( grid.mask, [1 1 nz] );
		mask = permute ( mask, [3 1 2] );
	end
	data = data .* mask;
end

transect_length = length ( vslice_x_pt );

%
% Geographic coordinates.
% find the cumulative distance along the transect in
% geographic coordinates.
switch gdata.rslice_grid.coordsys
case 'geographic'
	r = geodist_wr ( vslice_x_pt, vslice_y_pt );
	%r = m_lldist ( vslice_x_pt, vslice_y_pt );
	%r = distance_wr ( vslice_x_pt, vslice_y_pt );
case 'xy'
	r = sqrt ( diff(vslice_x_pt).^2 + diff(vslice_y_pt).^2 );
end


x = [0; cumsum(r)];

%
% Convert to kilometers
x = x/1000;
x = repmat ( x', nz, 1 );
transect_dist = x(1,:);







%
% Now interpolate the data.
switch ( gdata.interpolation_method )
case 'nearest_neighbor'

	for j = 1:length(vslice_x_pt)
	
	
		%
		% Assume the projection is ok to do this.
		d = sqrt ( (grid.x - vslice_x_pt(j)).^2 + (grid.y - vslice_y_pt(j)).^2 );
		ind = find(d == min(d(:)) );
		closest_point_index(j,1) = ind(1);
	
		%
		% we also have to find the closest rho points.  This may be the same.
		xd = (gdata.rslice_grid.x.data - vslice_x_pt(j));
		yd = (gdata.rslice_grid.y.data - vslice_y_pt(j));
		d = sqrt ( xd.^2 + yd.^2 );
		ind = find(d == min(d(:)) );
		closest_point_rho_index(j,1) = ind(1);
	
	
	
	
	end
	


	%
	% Assign the points
	[r,c] = indexit ( grid.x, closest_point_index );
	n = length(r);
	if strcmp ( z_type, '2D' )

		for j = 1:n
			return_data(j,1) = data(r(j),c(j));
		end

		z = [];
		return

	else % 3D


		return_data = zeros ( nz, transect_length );
		for j = 1:n
			tmp = squeeze(data(:,r(j),c(j)));
			return_data(:,j) = data(:,r(j),c(j));
		end
		return_data = flipud ( return_data );

	end


	%
	% Now create the y array (for display purposes).  In real physical
	% terms, it is actually z, or depth.
	[z,sc,Cs]=scoord_xy( gdata.rslice_grid.roms_grid.h, ...
	                     gdata.rslice_grid.roms_grid.theta_s, ...
			     gdata.rslice_grid.roms_grid.theta_b, ...
			     gdata.rslice_grid.roms_grid.Tcline, ...
			     length(gdata.rslice_grid.roms_grid.sc_r), ...
			     gdata.z_type, ...  
			     closest_point_rho_index );


case 'griddata_lite'
	
	xg = grid.x;
	yg = grid.y;
	%xg = gdata.rslice_grid.x.data;
	%yg = gdata.rslice_grid.y.data;
	z_x = gdata.rslice_grid.x.data;
	z_y = gdata.rslice_grid.y.data;
	z_r = gdata.rslice_grid.roms_grid.z_r;
	z_w = gdata.rslice_grid.roms_grid.z_w;

	%
	% This method relies upon a global variable called TRI.  We need
	% to clear it each time this is called.
	clear global TRI;
	global TRI;


	%
	% Compute the triangulation exactly once, the first time thru.
	if nz == 1
		%
		% Just a 2D field.
		ti = griddata_lite ( xg, yg, data, vslice_x_pt, vslice_y_pt );
		return_data = ti;
		z = zeros(size(ti));
	else
		%
		% 3D field.  Have to interpolate depths as well.
		clear global TRI;
		global TRI;
		ti(1,:) = griddata_lite ( xg, yg, squeeze(data(1,:,:)), vslice_x_pt, vslice_y_pt );
		for zi = 2:nz
			ti(zi,:) = griddata_lite ( xg, yg, squeeze(data(zi,:,:)), vslice_x_pt, vslice_y_pt );
		end



		%
		% Is it an s_r variable or an s_w variable?
		switch ( gdata.z_type )
		case 's_rho'
			zg = z_r;
		case 's_w'
			zg = z_w;
		otherwise
			error ( '%s:  unknown z type %s\n', mfilename, gdata.z_type );
		end



		%
		% Grid the proper depths.
		clear global TRI;
		global TRI;
		z_rw_i(1,:) = griddata_lite ( z_x, z_y, squeeze(zg(1,:,:)), vslice_x_pt, vslice_y_pt );
		for zi = 2:nz
			z_rw_i(zi,:) = griddata_lite ( z_x, z_y, squeeze(zg(zi,:,:)), vslice_x_pt, vslice_y_pt );
		end

		return_data = ti;
		z = flipud(z_rw_i);
	end


otherwise
	msg = sprintf ( '%s:  unhandled interpolation method.\n' );
	rslice_error ( msg );
end



%
% The last thing to do is to NaN out any of the data points that fell outside the
% grid.
n = length(return_data);
if strcmp ( z_type, '2D' )

	for j = 1:n
		if ~vslice_pt_inside_grid(j)
			return_data(j,1) = NaN;
		end
	end

	z = [];

else

	for j = 1:n
		if ~vslice_pt_inside_grid(j)
			return_data(:,j) = NaN;
		end
	end
	return_data = flipud ( return_data );

end




